\documentclass[12pt, a4paper]{article}
\usepackage{geometry}
\geometry{
	a4paper,
	left=12.7 mm,
	right=12.7 mm,
	top=12.7 mm,
	bottom=12.7 mm,
}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{physics}
\graphicspath{ {./} }
\usepackage{ctex}

\newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\formula}[1]{\text{式} \ref{#1} }
\newcommand{\matrx}[1]
{
	\ensuremath
	{
		\left (
		\begin{matrix}
			#1
		\end{matrix}
		\right)
	}
}

\begin{document}
	\section{Laplacian方程}
	Laplacian方程具有如下形式：
	\begin{equation}
		u = u(x,y,z) \qquad \laplacian u = 0
	\end{equation}
	其中$\laplacian = \pdv[2]{u}{x} + \pdv[2]{u}{y} + \pdv[2]{u}{z} $。
	由于Laplacian方程不显含时，因此Laplacian方程的解是“静态的”。
	
	\subsection{一维Laplacian方程}
	一维Laplacian方程相当于一个二阶ODE:
	\begin{equation}
		u = u(x) \qquad
		\begin{cases}
			\pdv[2]{u}{x} &= 0 \\
			u(x_0) & = a \qquad \text{边界条件1}\\
			u(x_n) & = b \qquad \text{边界条件2}\\
			x & \in (x_0, x_n) \\
		\end{cases}
	\end{equation}
	代入二阶导的离散形式，得到一个线性方程组：
	\begin{equation}
		\frac{u_{i-1} - 2 u_{i} + u_{i+1} }{(\Delta x)^2}= 0 \qquad i = 1,2,\dots,n-1
		\Rightarrow
		\begin{cases}
			u_0 &= a\\
			u_0 - 2u_1 +u_2 &= 0 \\
			u_1 - 2u_2 +u_3 &= 0 \\
			\dots \\
			u_{n-2} - 2u_{n-1} +u_{n} &= 0 \\
			u_n &= b \\
		\end{cases}
	\end{equation}
	尽管形式简单，但Laplacian方程与先前遇到的其他ODE不大相同：
	我们只知道$u$在两侧端点的值（而不是一侧端点及其导数值），这使我们无法直接迭代求解。
	比如，为了求解$u_2$，我们需要知道$u_3$，而$u_3$也是未知的。
	不过，眼尖的我们不难看出，这相当于这样一个矩阵方程：
	\begin{equation}
		\begin{bmatrix}
			1 & 0 & 0 & 0 & \cdots & 0 \\
			1 & -2 & 1 & 0 & \ddots & 0 \\
			0 & 1 & -2 & 1 & \ddots & 0 \\
			\vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\
			0 & \cdots & 0 & 1 & -2 & 1 \\
			0 & \cdots & 0 & 0 & 0 & 1
		\end{bmatrix}_{(n+1)\times (n+1)}
		\begin{bmatrix}
			u_0 \\ 
			u_1 \\ 
			u_2 \\
			\vdots \\ 
			u_{n-1} \\
			u_n
		\end{bmatrix}_{(n+1)\times 1}
		=
		\begin{bmatrix}
			a \\ 
			0 \\ 
			0 \\
			\cdots \\ 
			0 \\
			b
		\end{bmatrix}_{(n+1)\times 1}
		\Rightarrow
		L \bvec u = \bvec b
	\end{equation}
	其中$\bvec u$代表我们的场，$\bvec b$代表边界和内部“加载”项（此处为零），而矩阵$L$相当于二阶求导算符的离散形式：
	\begin{equation}
		L = \begin{bmatrix}
			1 & 0 & 0 & 0 & \cdots & 0 \\
			1 & -2 & 1 & 0 & \ddots & 0 \\
			0 & 1 & -2 & 1 & \ddots & 0 \\
			\vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\
			0 & \cdots & 0 & 1 & -2 & 1 \\
			0 & \cdots & 0 & 0 & 0 & 1
		\end{bmatrix}_{((n+1)\times(n+1) )}
	\end{equation}
	只需求解这样一个方程组，我们就能得到$\bvec u$。

	上述方程组更标准的做法是，将$u_0$和$u_n$排除出方程组并单独处理边界：
		\begin{equation}
		\begin{bmatrix}
			 -2 & 1 & 0 & \ddots & 0 \\
			 1 & -2 & 1 & \ddots & 0 \\
			 \ddots & \ddots & \ddots & \ddots \\
			 \cdots 0 & & 1 & -2 & 1  \\
			 \cdots 0 & & 0 & 1 & -2  \\
		\end{bmatrix}_{(n-1)\times (n-1)}
		\begin{bmatrix}
			u_1 \\ 
			u_2 \\
			\vdots \\ 
			u_{n-2} \\
			u_{n-1} \\
		\end{bmatrix}_{(n-1)\times 1}
		=
		\begin{bmatrix}
			0 \\ 
			0 \\ 
			\cdots \\ 
			0 \\
			0
		\end{bmatrix}_{(n-1)\times 1}
		+
		\begin{bmatrix}
			a \\ 
			0 \\ 
			\cdots \\ 
			0 \\
			b
		\end{bmatrix}_{(n-1)\times 1}
	\end{equation}
	这让$L$的形式更为对称。据我所知，两种写法是等价的。
	\newpage
	
	\subsection{Laplacian 方程}
    \footnote{参考：Sauer《数值分析》}
	二维 Laplacian方程具有如下形式：
	\begin{equation} \label{eq_1}
	u = u(x, y), 
	\left \{
	\begin{aligned} 
		\pdv[2]{u}{x} + \pdv[2]{u}{y} =0\\
		u(x_1,y) &= f1(y) \qquad \text{边界条件1}\\
		u(x_n,y) &= f2(y) \qquad \text{边界条件2}\\
		u(x,y_1) &= f3(x) \qquad \text{边界条件3}\\
		u(x,y_n) &= f4(x) \qquad \text{边界条件4}\\
		x & \in (x_1, x_n)\\
		y & \in (y_1, y_n)\\
	\end{aligned}
	\right.
	\end{equation}
	依旧离散化该方程（假定$\Delta x = \Delta y$）：
	\begin{equation}
	\begin{aligned}
		& \laplacian{u} = 0 \\
		\Rightarrow & \frac{u(x+\Delta x, y) - 2u(x, y) + u(x - \Delta x, y)}{(\Delta x)^2} + \frac{u(x, y+\Delta y) - 2u(x, y) + u(x, y-\Delta y)}{(\Delta y)^2}  = 0\\
		\Rightarrow & - 4 u(x,y) + u(x+\Delta x, y) + u(x - \Delta x, y)+ u(x, y+\Delta y) + u(x, y-\Delta y) = 0  \\
	\end{aligned}
	\end{equation}
	亦即对于$u$的内部网格点
	\begin{equation} 
		-4 u_{i, j} + u_{i+1,j} + u_{i-1,j}+ u_{i,j+1} + u_{i,j-1} = 0 \qquad i,j = 2,3,\dots,n-1
	\end{equation}
	$i,j$下标代表空间的离散步；为了更好地配合octave/matlab的规范，我们假定$u$的下标从$1$开始，而不是从$0$开始。
	和一维情形类似，我们依然无法使用迭代法解决这个问题。

	\newpage
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.3 \linewidth]{Repmat1}
		\caption{把$u$“拉直”为一个头尾相连的向量$\_u$ （为了方便描述，示意图取$n=3$）}
		\label{fig:repmat1}
	\end{figure}
	
	为了沿用矩阵方程组的思路，我们先把$u$“拉直”（使用Repmat 函数等）为一个头尾相连的向量$\_u$\footnote{此处使用Octave的规范。在不同软件、书本中，具体顺序可能不同}。
	如果$u$是一个$n\times n$的矩阵，那么 $\_u$是一个$n^2 \times 1$的向量，
	同时$u$的第$i$行第$j$个元素对应$\_u$的第$i+(j-1) n$行的元素，即$u_{i,j}$对应 $\_u_{i+(j-1) n}$等，
	以下设 $idx = i+(j-1)n$。
	如果使用Matlab，可以使用$sub2ind$等函数自动将$u$的$i,j$换算为$\_u$的$idx$，以增加可拓展性。
	$$
	u \rightarrow \_u
	$$
	随后，我们将 \formula{eq_1} 改写为关于 $\_u$的矩阵乘法形式：
	\begin{equation}
		A \cdot \_u = b
	\end{equation}
	显然，$A$是$n^2\times n^2$的矩阵，且$b$是一个$n^2 \times 1$的向量。
	我们讨论$A$的构造方法：
	\begin{itemize}
		\item 	对于处于$u$边界的点（$i==1, i==n, j==1, j==n$），
		由于$u_{i,j}$完全由边界确定，
		\begin{equation}
			\text{边界：} \quad u_{i,j} = f(x_i, y_i)
		\end{equation}
		因此直接取$A_{idx,idx} = 1, b_{idx} = f(x_i, y_j)$，$A$同行其余项为$0$：
		
		\begin{figure}[h]
			\centering
			\includegraphics[width=0.4 \linewidth]{Repmat2}
			\caption{边界点}
			\label{fig:repmat2}
		\end{figure}
		
		其中$f$是已知边界条件。
		\newpage
		
		\item 	
		对于处于$u$内部的点，
		$u_{i,j}$由迭代格式确定，
		\begin{equation}
			\text{内部：} -4 u_{i, j} + u_{i+1,j} + u_{i-1,j}+ u_{i,j+1} + u_{i,j-1} = 0
		\end{equation}
		因此取
		$A_{idx, idx} = -4,
		A_{idx, idx-1} =  
		A_{idx, idx+1} = 
		A_{idx, idx-n} = 
		A_{idx, idx+n} = 1,
		b_{idx} = 0
		$
		，$A$同行其余项为$0$：
		
		\begin{figure}[h]
			\centering
			\includegraphics[width=0.6 \linewidth]{Repmat3}
			\caption{内部点}
			\label{fig:repmat3}
		\end{figure}
	\end{itemize}
	可见，此时$A$相当于二维下$\laplacian$算符的离散形式。
	至此，我们原则上可以求解$\_u$：
	\begin{equation}
		\_u = A^{-1} b
	\end{equation}
	但由于$A$很大，并且大部分项是$0$，因此直接求逆会消耗大量计算资源和时间。
	聪明的做法是，使用稀疏矩阵存储$A$，并且使用相应的稀疏求解器，
	例如稀疏矩阵的GMRES(Octave和Matlab都自带该函数)：
	\begin{equation}
	\_u = \mathrm{GMRES}(A, b)
	\end{equation}
	然后再将$\_u$拓展为原先的矩阵
	$$
	\_u \rightarrow u
	$$
	至此，我们完成了Laplacian方程的求解。
	
\end{document}